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IN THE PROTOPLANETARY NEBULA 

Jeffrey N. Cuzzi 1 ' 4 , Robert C. Hogan 2 , and Karim Shariff 3 

Astrophys. J., Accepted, April 11 2008 

ABSTRACT 

We outline a scenario which traces a direct path from freely-floating nebula particles to the first 
10-100km-sized bodies in the terrestrial planet region, producing planetesimals which have properties 
matching those of primitive meteorite parent bodies. We call this primary accretion. The scenario 
draws on elements of previous work, and introduces a new critical threshold for planetesimal for- 
mation. We presume the nebula to be weakly turbulent, which leads to dense concentrations of 
aerodynamically size-sorted particles having properties like those observed in chondrites. The frac- 
tional volume of the nebula occupied by these dense zones or clumps obeys a probability distribution 
as a function of their density, and the densest concentrations have particle mass density 100 times 
that of the gas. However, even these densest clumps are prevented by gas pressure from undergo- 
ing gravitational instability in the traditional sense (on a dynamical timescale). While in this state 
of arrested development, they are susceptible to disruption by the ram pressure of the differentially 
orbiting nebula gas. However, self-gravity can preserve sufficiently large and dense clumps from ram 
pressure disruption, allowing their entrained particles to sediment gently but inexorably towards their 
centers, producing 10-100 km "sandpile" planetesimals. Localized radial pressure fluctuations in the 
nebula, and interactions between differentially moving dense clumps, will also play a role that must 
be allowed for in future studies. The scenario is readily extended from meteorite parent bodies to 
primary accretion throughout the solar system. 

Subject headings: solar systemfformation; accretion disks; minor planets: asteroids; turbulence; insta- 
bilities 
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1. INTRODUCTION 

There is no currently accepted scenario for the forma- 
tion of the parent bodies of primitive meteorites which 
accounts for the most obvious of their properties. These 
properties (reviewed by Scott and Krot 2005 and dis- 
cussed in more detail in section 2.1) include (a) dom- 
inance by aerodynamically well-sorted mineral particles 
of sub-mm size; (b) class-to-class variation in well-defined 
physical, chemical, and isotopic properties; (c) a spread 
of 1 Myr or so between the formation times of the old- 
est and youngest objects found in the same meteorite; 
(d) a spread of 1-3 Myr in radiometric ages of different 
meteorite types; and (e) a dearth of melted asteroids, 
with model results for even some melted asteroids which 
imply Myr delays in formation relative to ancient miner- 
als. In recent years, meteoritic evidence has appeared for 
some very early-formed planetesimals, which represent a 
minority of both meteorites and asteroids. This implies 
that primary accretion started early and continued for 
several million years - thus, it was fairly inefficient and 
did not run quickly to completion. 

Most current models for this "primary accretion" stage 
(reviewed by Cuzzi and Weidenschilling 2006 and Do- 
minik et al 2007) can be classified as either (a) incre- 
mental growth, where large particles sweep up smaller 
ones by inelastic collisions involving porous surfaces, and 
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growth proceeds hierarchically; or (b) instability, where 
physical sticking is irrelevant and collective effects drive 
collapse to km-sizes or larger on very short timescales. 
Those who favor instability models, most of which rely 
on gravity and occur in a particle-rich nebula midplane, 
are concerned by the poorly understood sticking of min- 
eral particle aggregates and the apparent difficulty of 
growing beyond meter size due to rapid inward migration 
and collisional disruption. Those who favor incremental 
growth have noted that midplane instability models are 
precluded by even very weak global turbulence, and that, 
in the dense midplane layers that form when turbulence 
is absent, incremental accretion is at low relative veloc- 
ity and the meter-size barrier is not a problem. However, 
the most sophisticated models of incremental accretion in 
nonturbulent nebulae find it to be so efficient that large 
planetesimals grow in only 10 4 — 10 5 years throughout 
the asteroid belt region (Weidenschilling 2000), a short 
time which is difficult to reconcile with constraints (a-e) 
above. A third class of scenario suggests that a com- 
plex interplay between several nonlinear processes - tur- 
bulence, pressure gradients, and gravity - may concen- 
trate appropriately sized particles and lead to planetes- 
imal growth (Cuzzi et al 2001, 2005, 2007; Johansen et 
al 2006, 2007). Scenarios which involve preferentially 
meter-sized objects encounter concerns about whether 
meter-sized objects can survive their own high-velocity 
mutual collisions in this sort of turbulent environment 
(cf. Sirono 2000, Langkowski et al 2007, Ormel et al 
2007). However, in principle, such boulder-concentration 
scenarios can proceed independently in the same envi- 
ronment as discussed in this paper, where we focus on 
mm-size particles. 
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Most of the above scenarios provide no natural ex- 
planation for observed meteorite properties (a) and (b) 
above. In particular, the evidence for the H chondrite 
class (and perhaps all the ordinary chondrites) suggests 
that entire asteroids of 10-100km diameter formed di- 
rectly from a physically chemically and isotopically ho- 
mogeneous mix of dust-rimmed particles of similar size 
(section 2.1). In this paper, we outline a possible path by 
which entire batches of mm-size, aerodynamically sorted 
particles might proceed directly in turbulence (even if 
sporadically) to planetesimals having the properties out- 
lined above. We find that sufficiently large and dense 
chimps of mm-size particles can form by turbulent con- 
centration such that, even if classical gravitational in- 
stability can't operate, their self-gravity may still allow 
them to survive disruptive forces and slowly sediment 
into a "sandpile" planetesimal. This simple analytical 
model is backed up by some numerical simulations that 
support the basic idea. In sections 2.2 and 2.3, we review 
the most relevant physics that determines the properties 
of dense, particle-rich zones or clumps in turbulent neb- 
ulae. In section 3, we address the fate of these dense 
clumps using analytical and numerical models of their 
evolution. We derive the combination of size and mass 
density a clump must have to evolve into a "sandpile" 
having some degree of internal strength. In subsequent 
stages not modeled here, we imagine that collisions and 
thermal sintering transform these sandpiles into the co- 
hesive rocky parent bodies we see today. However, in the 
outer solar system, lower energy collisions and weaker 
thermal processing might well allow planetesimals to re- 
tain their initial low-strength states, as seen for some 
cometary objects {eg., Asphaug and Benz 1996). 

2. BACKGROUND 

2.1. Meteoritics background and evidence for inefficient 
accretion in a turbulent environment 

One can extract a number of clues from primitive (un- 
melted) meteorites regarding the primary accretion pro- 
cess by which their parent bodies first formed (Scott 
and Krot 2005, Taylor 2005). For the best evidence, 
one must look back through extensive subsequent evo- 
lutionary stages. Even unmelted bodies in the 100km 
size range have incurred extensive collisional evolution 
(Bischoff et al. 2006), producing compaction, fragmen- 
tation, and physical grinding and mixing on and beneath 
their surfaces, which may obscure the record of primary 
accretion. Model studies (e.g., Petit et al 2001, Kenyon 
and Bromley 2004, 2006; Bottke et al 2005; Chambers 
2004,2006; Weidenschilling and Cuzzi 2006) suggest that 
the collisional stage occurred after dispersal of the nebula 
gas allowed the orbital eccentricities of primitive bodies 
to grow. We are concerned with an earlier stage, when 
the still-abundant nebula gas led to a more benign en- 
vironment with fewer and gentler collisions. The direct 
products of primary accretion might be most clearly visi- 
ble in the rare "primary texture" seen in some CM (Met- 
zler et al 1992) and CO (Brearley 1993) chondrites. In 
these objects, or more specifically in unbroken fragments 
within them, the texture consists of similarly-sized, dust- 
rimmed particles packed next to each other as if gently 
brought together and compressed, with no evidence for 
local fracturing or grinding. 

Even after collisional effects associated with subse- 



quent stages of growth have blurred this signal, perhaps 
even mixing material from different parent bodies, ev- 
idence remains in the bulk properties of all chondrite 
classes. Chondrite classes are defined by their distinctive 
mineral, chemical, and isotopic properties (e.g. Gross- 
man et al 1988, Scott and Krot 2005, Weisberg et al 
2006). Large samples of material with a quite well de- 
fined nature were accumulated at one place and/or time, 
and material of a quite different, but equally well-defined, 
nature was accumulated at another place and/or time 
into a different parent body. Amongst the most obvious 
aspects of these class properties is a dominance within 
chondrites of sub-mm size mineral particles (generically 
"chondrules" ; Grossman 1988, Jones et al 2000, 2005; 
Connolly et al 2006) which are aerodynamically well- 
sorted. Aerodynamic sorting of chondrite components 
was first emphasized by Dodd (1976), Hughes (1978, 
1980), and Skinner and Leenhouts (1993), and has since 
been discussed by a variety of authors (see Cuzzi and 
Weidenschilling 2006 for a summary of the observations 
and arguments supporting aerodynamics). Chondrules 
have a size which varies from class to class, but is distinc- 
tive and narrowly defined within a given class or chon- 
drite. 

The evidence suggests that chondrules do not com- 
prise merely a thin surface layer swept up by a large 
object (Scott 2006). In the case which we suggest as an 
archetype, the H-type ordinary chondrite parent body 
is widely believed to be a 100 km radius object (per- 
haps the asteroid Hebe), initially composed entirely of 
a physically, chemically, and isotopically homogeneous 
mix of chondrules and associated material, which was 
thermally metamorphosed by accreted 26 Al to a degree 
which varied with depth, into an onion-shell structure, 
and subsequently broken up in several stages (Taylor et 
al 1987, Trieloff et al 2003, McSween et al 2002, Grimm 
et al 2005). Bottke et al (2005) conclude from models 
of collisional evolution in the primordial (massive) and 
the current (depleted) asteroid belt, that the primordial 
asteroid mass distribution was dominated by objects hav- 
ing diameter of around 100km, rather than having a pow- 
erlaw size distribution somewhat like that of the current 
population. 

We see the essential challenge as understanding how 
such a large object can be assembled from such a ho- 
mogeneous mixture of mm-size constituents, while other 
objects (arguably the parents of the L and LL type ordi- 
nary chondrites as well, and logically then the parents of 
the chondrites of all classes) are assembled from distinct, 
yet comparably homogeneous, ensembles of qualitatively 
similar particles. Moreover this assembly, or primary ac- 
cretion, phase of nebula evolution must persist for a du- 
ration of lMyr or more between the formation times of 
the oldest and youngest meteorites (Wadhwa et al 2007) 
and even of objects found in the same meteorite. The > 
1 Myr age spread between ancient refractory inclusions 
(CAIs) and chondrules is well known (Russell et al 2006), 
but there may even be a comparably extended age spread 
between different chondrules in a given meteorite (Kita 
et al 2005). 

In recent years, radiometric dating of achondrites 
(melted meteorites) has revealed some early-formed (and 
early melted) planetesimals (Kleine et al 2005, 2006). 
One expects early formation and melting to go together, 
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because objects larger than only 10km or so, which ac- 
creted at the same time that refractory inclusions formed 
(4.567 Byr ago), would have incorporated their full com- 
plement of 26 Al and would have melted (LaTourette and 
Wasserburg 1998; Woolum and Cassen 1999, Hevey and 
Sanders 2006). Thus, early accretion implies melting, 
and lack of melting requires late accretion. Model re- 
sults imply that, to avoid complete melting, 100km bod- 
ies must accrete only after a typically 1.5-2 Myr delay 
relative to ancient CAI minerals (reviewed by McSween 
et al 2002 and Ghosh et al 2006). While unknown se- 
lection and sampling effects may influence the limited 
meteorite data record, spectral analysis of the entire as- 
teroid database also suggests that melted asteroids are 
rare. That is, while debate continues as to whether small 
amounts of melt might be found on the surfaces of many 
asteroids, which might be caused by impacts (e.g. Gaffey 
et al 1993, 2002), only Vesta and a handful of other ex- 
amples of fully melted basaltic objects have been found in 
spite of extensive searches (Binzel et al 2002; Moskovitz 
et al 2007). 

Put together, this body of information suggests that 
primary accretion started early but continued for several 
million years - thus, it was inefficient, or at least spo- 
radic. By comparison, midplanc incremental accretion 
in a nonturbulent nebula accretes numerous Ceres size 
bodies and tens of thousands of 10-100km size objects in 
only 10 5 years (Weidenschilling 2000); thus, it is highly 
efficient. 

Models suggest that the density, temperature, and 
composition of the nebula evolved significantly over sev- 
eral million years (Bell et al 1995, D'Alessio et al. 2005, 
Ciesla and Cuzzi 2006). Thus one expects the properties 
of forming planctcsimals to change with time even though 
the chemical, isotopic, and mineralogical properties of 
the mixture of solids from which they formed might have 
been fairly uniform in space at any given time (see Cuzzi 
et al 2003 and Ciesla 2007 for more detailed discussions) . 
Reality might have been even more complicated - there is 
some evidence that chondrites of very distinct chemical 
and isotopic properties might have formed at roughly the 
same time (Kita et al 2005), suggesting a combination of 
spatial and temporal gradients. 

In this paper we present the overview of a scenario by 
which accumulation of a well-defined particle mix into a 
large planetesimal - making it a "snapshot" of the local 
nebula mix - might occur sporadically or inefficiently in 
weak turbulence, over an extended period of time and 
throughout the terrestrial planet region. The challenge 
for the future is to show that the accretion processes 
discussed here can quantitatively account for the mass 
needed to build planetesimals of the mass needed to cre- 
ate planetary embryos and planets during the age of the 
nebula (Chambers 2004, 2006; Cuzzi et al 2007). 

2.2. Turbulence 

A variety of astronomical and meteoritical evidence 
supports weak, but widespread and sustained, nebula 
turbulence. The observed abundance of small particles 
at high altitudes above the midplanc in many visible 
protoplanctary nebulae for millions of years is most eas- 
ily explained by ongoing turbulence - both to regenerate 
the small grains in collisions and to redistribute them to 
the observed altitudes (Dullemond and Dominik 2005). 



The survival of ancient, high-temperature mineral inclu- 
sions (CAIs) in chondrites, in spite of their tendency 
to drift inwards towards the sun, can also be explained 
by turbulent diffusion (Cuzzi et al 2003, 2005). More- 
over, the recent discovery by STARDUST of both high- 
and moderate-temperature crystalline silicates in comet 
Wild2 (Brownlee et al 2007, Zolensky et al 2007) can be 
easily understood in this way (Ciesla 2007). 

It remains a subject of active theoretical debate as to 
how, or even whether, the nebula can maintain itself in 
a turbulent state. Ongoing dissipation by molecular vis- 
cosity v m on the smallest, or Kolmogorov, lengthscale 
requires a production mechanism to continually generate 
turbulent kinetic energy. Although the very existence of 
ongoing disk accretion into the star, with the ensuing 
gravitational energy release, provides an ample energy 
source, the vehicle for transferring that energy into tur- 
bulent gas motions is undetermined. Baroclinic instabil- 
ities require substantial opacity at thermal wavelengths, 
but might be suitable to generate turbulence during the 
fairly opaque chondrule-CAI epoch in the inner solar 
system (Klahr and Bodcnhcimcr 2003). The widely ac- 
cepted magnetorotational instability (MRI) (Stone et al 
2000) may or may not be precluded at the high gas den- 
sities of the inner solar system (Turner et al 2007). The 
role of pure hydrodynamics operating on radial shear 
remains unclear, because Kcplerian disks are nominally 
stable to linear (small) perturbations (Balbus et al 1996, 
Stone et al 2000). However, a number of current theoret- 
ical studies suggest turbulence might indeed be present 
(eg. Arldt and Urpin 2004, Busse 2004, Umurhan and 
Regev 2004, Dubrulle et al 2005, Afshordi et al 2005, 
Mukhopadhyay et al 2005, Mukhopadhyay et al 2006). 
Much of this recent work focusses on nonlinear (finite am- 
plitude) instabilities occurring at the very high Reynolds 
numbers characterizing the nebula, which are well be- 
yond the capability of current numerical and even labo- 
ratory models. 

In this paper, we simply assume the presence of weak 
nebula turbulence. The intensity of turbulence may 
be characterized by a nondimensional parameter a — 
(vt/c) 2 , where Vt is the velocity of the most energetic 
eddies and c is the sound speed. Values of a in the 
range 10~ 4±1 are consistent with observed nebula life- 
times, and can be sustained by only a few percent of the 
gravitational energy released as the nebula disk flows into 
the sun. One then estimates the most energetic (largest) 
eddy size as L ~ Ha 1 / 2 , where H is the nebula verti- 
cal scale height, or disk thickness. A cascade to smaller 
scales ensues, where the smallest (Kolmogorov) eddy 
scale is r\ = LRc~ 3 / 4 ~ 1 km and Re = acH/v m <~ 10 7 
for a <~ 10~ 4 is the nebula Reynolds number (Cuzzi et 
al 2001). 

2.3. Turbulent concentration 

It has been found both numerically and experimen- 
tally that particles of a well defined size and density 
are concentrated into very dense zones in realistic 3D 
turbulence (Eaton and Fessler 1994). The optimally 
concentrated particle has a gas drag stopping time t s 
which is equal to the overturn time of the Kolmogorov 
scale eddies. In the nebula, t s = rp s /cp g where r and 
p s are the particle radius and material density, and c 
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Fig. 1. — Cumulative probability distribution, or volume frac- 
tion F p , of particles having local concentration factor larger than 
some value of C (from Cuzzi ct al 2001; see also section 2.3). The 
curves with symbols represent comparison of the predictive theory 
described in Cuzzi ct al (2001) with 3D numerical simulations at 
moderate Reynolds numbers; the curves without symbols are ex- 
tension of the same theory to nebula conditions denoted by their 
values of a. Concentration factors C ~ 100 of chondrule precur- 
sors (leading to <3? = p p /p g ~ 1) are fairly common for a range of 
nebula turbulent intensities (denoted by a), and help explain lack 
of chondrule isotopic fractionation (Cuzzi and Alexander 2006). In 
this paper we focus on the much less common, but much denser, 
regions towards the bottom right side of the plot (C ~ 10 4 or 
<E> ~ 100). More general PDFs which separate out the role of 
vorticity and incorporate the effect of particle mass loading are 
discussed in section 2.3 and Appendix A. 

and pg are the gas sound speed and density. The Kol- 
mogorov eddy frequency cu n scales with the orbital fre- 

quency fl as to v = f2Re ' , and particles are most ef- 
fectively concentrated when t s u v ~ 1 (Eaton and Fessler 
1994). Chondrule-size silicate particles are concentrated 
for canonical nebula properties if a ~ 10~ 5 — 10~ 3 ; more- 
over the shape of the concentrated particle size distribu- 
tion is parameter- independent, and agrees very well with 
that of chondrules (Paque and Cuzzi 1997, Cuzzi et al 
2001). We define the concentration factor C = p p /p~p~, 
where we denote the local volume mass density of parti- 
cles as p p and p~p~ is its nebula average value, and the par- 
ticle mass loading $ = p P /p g — C~pp~/p g . For instance, 
figure 1 (Cuzzi et al 2001) shows the probability distri- 
bution function (PDF) for C determined from numerical 
simulations at moderate Reynolds number (curves with 
symbols) and predicted for plausible nebula Reynolds 
numbers (curves without symbols). Note that the vol- 
ume fraction having C > 100 is quite substantial, which 
helps explain some mineralogical and isotopic properties 
of chondrules (Cuzzi and Alexander 2006). 5 Not surpris- 
ingly, the volume fraction decreases for larger C. In this 
work we focus on the less abundant, but higher mass den- 

5 Here, <I> = p p /p g is different from the definition in Cuzzi and 
Alexander (2006), which is equivalent to p p /Ap g , where A ~ 10~ 2 
is the fractional abundance of solids to hydrogen gas by mass. Thus 
values of $ in Cuzzi and Alexander (2006) arc the equivalent of 
Pp/pg ~ 1 - a much more common occurrence than the situation 
we study here. 



sity, concentrations towards the lower-right hand side of 
the PDF. 

The turbulent concentration PDFs of figure 1 did not 
allow for the feedback of the concentrated particles on 
the gas. This "mass loading" would be expected to damp 
gas turbulent motions once the particle mass density sig- 
nificantly exceeds that of the gas, and lead to a satura- 
tion of concentration at some value. Recent work which 
includes the feedback effects, through drag, of particle 
density on the damping of turbulence, has shown that 
values of $ < 100 can indeed be achieved - although less 
commonly than in the unloaded models such as shown in 
figure 1. This cascade model of turbulent concentration 
(Hogan and Cuzzi 2007) is summarized in Appendix A. 
Based on these results, we adopt $ = 100 as an upper 
bound in all stability and evolution modeling. 

3. PRIMARY ACCRETION OF PLANETESIMALS 

The ability of gravity to overwhelm opposing forces is 
a well-known theme in astrophysics, dating back 80 years 
to early studies of star formation by Jeans. Over 30 years 
ago, gravitational instability (CI) of solids in a dense 
layer near the nebula midplane was proposed to lead di- 
rectly to solid planetesimals on a dynamical timescale 
(the orbit time) (Safronov 1969, Goldreich and Ward 
1973). Traditional CI thresholds require the dynami- 
cal (collapse) time of a dense region to — 7r(Gp p ) -1 / 2 , 
where G is the gravitational constant, to be shorter than 
both the transit time due to random velocity and the 
local shear (vorticity) timescale (Toomre 1964). These 
arguments imply that GI occurs when the local par- 
ticle density exceeds a few times the Roche density 
Pr = 3M /47ra 3 , where M Q is the solar mass and a 
is the distance from the sun (Safronov 1991, Cuzzi ct 
al 2001). However, this scenario is frustrated by self- 
generated midplane turbulence even if the nebula is not 
globally turbulent (e.g. Wcidenschilling 1980,1995; Cuzzi 
et al 1993,1994; Johansen et al 2007). That is, a dense 
layer of cm-size and larger particles orbits at a differ- 
ent velocity than the pressure-supported gas, leading to 
a vertical velocity shear which becomes turbulent and 
stirs the particle layer. A more recent flavor of GI re- 
quires the particles to be small enough that the particle- 
rich midplane acts as a single fluid, like fog, so that its 
vertical density stratification stabilizes it against self- 
generated turbulence (Sekiya 1998, Sekiya and Ishitsu 
2001, Youdin and Shu 2002, Sekiya and Takeda 2003, 
Youdin and Chiang 2004). Because the particles must be 
very small to satisfy this one-phase requirement (mm-size 
and smaller), their settling towards the midplane is frus- 
trated by even the faintest breath of nebula turbulence 
(a ~ 10~ 9 ). If the particles are prevented from settling 
into a sufficiently dense layer that their mass density can 
stratify the total density (eg, Dubrulle et al 1995, Cuzzi 
and Weidenschilling 2006), the conditions needed to trig- 
ger the instability are never achieved. 

In the following subsections, we study the evolution of 
a dense clump such as might form towards the lower right 
domain of figure 1, albeit in a simplified way. We treat 
a clump as an isolated object when in reality it exists in 
a context of surrounding material of varying density. In 
section 4, we discuss the implications of more realistic 
conditions. 
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Fig. 2. — Results of a ID compressible numerical model, show- 
ing how gas pressure precludes dynamical timescale gravitational 
instability (Sckiya 1983) as described in more detail in section 3.1. 
Here, Ma = (//c)(Gp 9 ) 1/2 = 1/$* is like a Mach number. Neb- 
ula values of c and p g lie in the top set of nearly linear curves, 
allowing only slow inward sedimentation on the timescale t se( j (for 
"£>() = 100). That is, under nebula conditions of gas density and 
temperature, and particle loading, the gas behaves as if it were 
incompressible and prevents dynamical timescale GI for particles 
with stopping time much less than the dynamical collapse time. 

3.1. Gas pressure precludes particle gravitational 
instability 

The role of the gas pressure in GI for nebula parti- 
cles has been almost universally ignored, even though 
its importance was pointed out decades ago by Sekiya 
(1983). We encountered it ourselves independently. To 
test our GI thresholds, we initially ran numerical models 
of static, but very dense, clumps (<E> = 1000), expecting 
them to collapse with gravitational free-fall or dynami- 
cal collapse times tG = Tr{Gp p )~ x / 2 = n(G^p g )^ 1 ^ 2 and 
velocities Vq w l(Gpp) 1 ^ 2 /2tt, where I is the clump di- 
ameter. We assumed a spherically symmetrical dense 
particle clump with density p p , embedded in gas with 
density p g , where $ = p P /p g - Instead, only very slow 
shrinkage ensued. We suspected that the collapse was 
being artificially blocked by our incompressible code. 

We then developed a fully compressible ID model 
which confirmed our incompressible calculations and 
clearly demonstrated a much higher, gas-pressure- 
dependent <E> threshold for GI in the limit when t s <C £g- 
Results from the model are shown in figure 2. 

Figure 2 shows results from this fully compressible 
model, illustrating how $ > c/(/(Gp 9 ) 1 / 2 ) is required 
for traditional (dynamical timescale) instability to occur. 
The initial condition was a Gaussian blob having den- 
sity distribution p p (r) — <f>op g cxp(— r 2 // 2 ). Slow shrink- 
age of the particles through the gas is seen in the sta- 
ble, "incompressible" regime where nebula parameters 
lie, well within the flat curves at the top of the plot with 
Ma - 10~ 4 for I ~ 10 4 km, c - 10 5 cm/s, and p g ~ 10~ 9 g 
cm~ 3 . The high value of Ma required for true dynamical 
timescale collapse (given the chosen value of $ = 100) 
would require smaller gas sound speed or larger clump 
size I and gas density p g , by orders of magnitude. 

The slow shrinkage in the stable regime results from 
particles sedimenting inwards towards their mutual cen- 



ter under their own self-gravity at their terminal velocity 
vt = gt s , where g = 2G$lp g is the local gravitational ac- 
celeration due to the clump's mass, and t s = rp s /cp g is 
the particle stopping time. Thus vt = 2G$lrp s /c <C Vq 
for t s <C to, and the "sedimentation" timescale is 

I 1 _c _L 

tsed ~2v~^~ 4G$(rp s /c) ~ 4G^rp s ~ 4G1>p g t s ' 

(1) 

roughly 30-300 orbit periods at 2.5 AU for 300p radius 
chondrules and $ = 1000 — 100; that is, much longer than 
to- It seems inappropriate to describe this slow ongoing 
evolution as an instability - certainly in midplane scenar- 
ios where it was slow vertical sedimentation by individual 
particles toward the midplane that produced the "un- 
stable" situation in the first place, and where only fur- 
ther slow sedimentation transpires. The situation is more 
akin to star formation mediated by ambipolar diffusion 
than by traditional gravitational instability. Gas pres- 
sure even precludes the Safronov-Goldreich-Ward GI in 
the form it was originally proposed, where cm-sized par- 
ticles were envisioned (although t se d would be faster for 
literally cm-sized particles, and more closely approaches 
to). 

The results of Sekiya (1983) have apparently been over- 
looked by all subsequent workers, who have invariably as- 
sumed a Toomre-Safronov-Goldreich-Ward type particle 
ensemble which is decoupled from the gas, and stabilized 
(on small scales) only by particle random velocities. This 
gas pressure constraint is, however, crucial for particles 
in the chondrule size regime where t s ~ an hour and 
to ~ a year, and must be applied to all GI models in- 
volving chondrule-sized particles. The good news is that 
for the interesting range of nebula gas density and par- 
ticle concentration values, the gas does actually behave 
as if it were incompressible (validating the use of an in- 
compressible code in numerical models such as shown in 
section 3.4). 

The physics is easily understood in the one-phase 
regime (t s <C to)- The particles feel no direct pres- 
sure force, but are firmly trapped to the gas which does. 
The (inward) gravitational force is dominated by parti- 
cles for $ ^> 1. The particles drag and compress the 
gas with them as they begin to collapse under their self- 
gravity, producing a radial gas density gradient dp g /dl; 
this translates into a large outward pressure gradient 
c 2 dp g /dl which acts on the gas, and thereby also on the 
strongly coupled particles. 

The gravitational force per unit volume is fa — 
AGMp p /l 2 — 2G$ 2 p 2 l, where M is the clump mass and 
R its radius. After some incipient shrinkage and com- 
pression occurs, the outward pressure gradient force per 
unit volume becomes roughly fp = c 2 dp g /dl ~ c 2 p g /l, 
where c is the gas sound speed. Then for GI to be possi- 
ble, fc/fp = G$ 2 l 2 p g /c 2 > 1, giving the criterion for GI 
as $ > $*, where $* = c/Z(Gp g ) 1//2 . <&* is like an inverse 
Mach number Ma" 1 , where Ma = (1/ c)(Gp g ) 1/2 . An al- 
ternate (perturbation) approach compares the incremen- 
tal changes Sf p and 5 fa associated with a small gravi- 
tationally induced shrinkage characterized by 51 I. It 
is straightforward to show (for t s <C to) that Sf p /5fa ~ 
c 2 p g /GMQ ~ 10 5 under nebula conditions. Thus the 
clump is "stiffened" by pressure and even a small incre- 
mental shrinkage is self-limiting. Inspection of figure 
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2 shows small oscillations which suggest that a formal 
linear stability analysis might lead to an improved sta- 
bility criterion; this is, however, beyond the scope of the 
present paper. 

For I ~ a few 10 3 km and p g ~ 10~ 9 g cm~ 3 (De- 
sch et al 2005), $* - 10 4 . This criterion for GI can 
be compared with the traditional criterion for (marginal, 
or maximal) instability which is widely used as the crite- 
rion for gravitational instability, roughly twice the Roche 
density p R = 3M Q /47ra 3 (Goldreich and Ward 1973, 
Safronov 1991, Cuzzi et al 1993). The ratio of the two 
thresholds is: 

<5>*p g p g 2c{Gp g ) 1/2 2 - 6 x 10 3 
2p R ~ 2p R M& ~ fWl (Z/10 3 km) ' 

for p g = 10~ 10 to 10~ 9 . For the midplane GIs which have 
been widely discussed in the past (reviewed by Cuzzi 
and Wcidcnschilling 2006), I ~ 10 2 km, so the traditional 
GI criterion falls short by a factor of more than 10 4 
(cf. Sekiya 1983)! This degree of mass concentration 
is unlikely under any plausible circumstances, especially 
for particles which are already well trapped to the gas. 
For the ubiquitous dense clumps we discuss here, which 
form at all elevations, I ~ 10 4 km. Thus the pressure- 
supported gas phase prevents the tightly coupled parti- 
cles from undergoing GI until the particle mass loading is 
hundreds of times larger than the traditional GI criterion 
(which is on the order of the Roche density p R ). 

Nevertheless, Sekiya (1983) also noted that when p p 
is in the range usually cited for GI (a few times p R , or 
$ ~ a few hundred), a 3D "incompressible mode" arises 
which, while not explicitly stated, we interpret as retain- 
ing an identifiable cohort of particles. An "incompress- 
ible mode" might resemble a blob of water oscillating 
in zero-gravity. Dense zones which result from turbu- 
lent concentration are candidates for such incompressible 
modes, even at mass loadings too low to induce actual GI. 
However, the fate of such slowly evolving entities in the 
presence of likely perturbations has never been explored. 
Below we discuss the least avoidable, most pervasive dis- 
ruptive perturbation which such clumps would encounter 
once they form, which we believe to be ram pressure dis- 
ruption due to systematic velocity differences between 
the dense zones and the gas, and then derive conditions 
under which these clumps can survive to become plan- 
etesimals. 

3.2. The fate of "incompressible" clumps 

Several kinds of perturbation by the enveloping nebula 
gas might disrupt a dense clump before the slow sedimen- 
tation of its constituent particles towards their mutual 
center (on the timescale t se d) can produce a moderately 
compact object. 

Perhaps the simplest to discuss and dismiss are tur- 
bulent pressure, velocity, and/or vorticity fluctuations 
encountered by the clump. Turbulent pressure fluctua- 
tions on the scale I have typical intensity p g v t (l) 2 , where 
Vt{l) is certainly smaller than the fluctuating velocity 
v t (L) = ca 1 ' 2 of the largest eddy (the dominant energy- 
containing eddy). The ratio of even the largest eddy 
pressure fluctuations to the steady ram pressure from 
the headwind /3Cla = (3vk (a is the semimajor axis) 
is (v t /j3v K ) 2 = ac 2 /(3 2 v 2 K . Since c/v K = H/a and 



P = (H/a) 2 , (v t /(3v K ) 2 < a/[3 < 1 unless a > 10~ 3 . 
Nebula evolutionary timcscalcs, and parameters leading 
to turbulent concentration of chondrule-sized particles 
in the asteroid belt region, imply 10~ 6 < a < 10~ 3 de- 
pending on gas density and other properties (Cuzzi et al 
2001). This suggests that turbulent velocity fluctuations 
on clump lengthscales play a small role. We have run 
simulations of clumps settling in gravity with and with- 
out enveloping turbulence, and the typical evolutions are 
qualitatively similar. 

We next address local average vorticities on the length- 
scale of a clump. It has long been known that local 
vorticity is a factor in gravitational instability ( Toomrc 
1964, Goldreich and Ward 1973). Because eddies much 
smaller than the largest scale of turbulence typically have 
larger vorticity than the Keplerian shear commonly ex- 
plored (eg. Toomre 1964), this is a concern in princi- 
ple. For example, in the inertial range, uj(1) ~ f2(L/Z) 2 / 3 
(Tennekes and Lumley 1972). Then for I ~ 10 4 km, 
L ~ Ha 1 / 2 and a ~ 10" 4 , uj(1) ~ 10 4 f2. However, turbu- 
lent concentration has the property that dense particle 
zones preferentially lie in zones of locally low vorticity 
compared to the average at their lcngthscale (Eaton and 
Fessler 1994). Statistical studies (Hogan and Cuzzi 2007) 
show that dense zones in high Re environments can form 
in regions with local vorticity 1-2 orders of magnitude 
smaller than the average value expected for that length- 
scale (see also figure 5). Below we sketch how such a 
constraint is derived and applied. 

A simplified requirement for gravitational binding of a 
clump of lcngthscale I and mass density p p , in the pres- 
ence of local rotation at angular velocity u>(l), is AGp v > 
lo 2 (1). Recalling that the local vorticity u(l) can vary 
by orders of magnitude from its average value (w(Z)), we 
normalize both sides of the above relationship between 
p p and u> by the average inertial range enstrophy (oj 2 (1)), 
and define the normalized enstrophy S = (u 2 (I) / (u 2 (I))) . 
This quantity provides the horizontal axis on figure 5. 
Then the simple expression above can be rewritten as 
$ > {uj 2 (l))S/4:Gpg. Inertial range relationships can be 
used to express (uj 2 (1)) = fl 2 (L/l) 4 ^ = n 2 (L/B V ) 4 / 3 
where r\ is the Kolmogorov scale and B is some scal- 
ing factor. Making use of Re = acH/v m = (L/rj) 4 / 3 
as discussed earlier, and substituting nebula parameters 
at 2.5 AU, we find that gravitational binding requires 
$ > 3 x 10 e aS. This constraint yields diagonal lines on 
plots of the sort shown in figure 5 (for an example sec 
Cuzzi et al 2007). Interesting regions of parameter space 
remain accessible to a degree that depends on quantita- 
tive modeling of the PDF P($, 5"). 

Obviously, more detailed numerical models of clumps 
in realistic turbulence are needed to assess these pertur- 
bations. We will assume, for the purposes of this study, 
that turbulent pressure (and associated vorticity) fluctu- 
ations are a minor influence and that regions which are 
stable to their own local vorticity exist (see Cuzzi et al 
2007). This leaves the dominant disruptive process, un- 
avoidably shared by turbulent and laminar nebulae, as 
ram pressure disruption due to systematic velocity dif- 
ferences between the clump and the gas, as described 
below. Incidentally, survival against ram pressure gives 
the horizontal stability boundaries in figure 3 of Cuzzi 
et al (2007). Emplacing these constraints on PDFs ob- 
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taincd from cascade models is somewhat involved, and a 
full quantitative description of the situation is deferred 
to a future paper. 

3.3. Ram pressure disruption and the Weber number 

Here, for simplicity, we envision the evolution of a sin- 
gle, isolated dense clump with $ 3> 1, merely as an ex- 
ample of one of the many dense clumps formed near the 
lower right hand corner of figure 1. The clump exerts 
a strong collective influence on the entrained gas, and 
the mixture experiences solar gravity as a unit. It at- 
tempts to move as an individual large object relative to 
the nebula gas, dragging the entrained gas along. For in- 
stance, a clump forming at some distance above the neb- 
ula midplane settles downward under the acceleration 
of the vertical component of solar gravity, whereas the 
surrounding gas remains under vertical hydrostatic bal- 
ance. Also, nebula gas generally orbits at a slower veloc- 
ity than the Keplerian orbit velocity vk obeyed by mas- 
sive particles. This is because of its generally outward 
radial pressure gradient force —{l/p g )dP/da = 2/3Cta 2 , 
where a is the distance from the sun, and (3 ~ 10 -3 ). 
Thus, particles experience a headwind with magnitude 
w g ~ [3Qa = /3vk ~ a few 10 3 cm/sec. These head- 
winds result in a ram pressure p g w 2 /2 which can disrupt 
a strengthlcss clump. 

For example, a dense drop of ink settling in a glass of 
water is quickly disrupted by the Rayleigh- Taylor insta- 
bility and mixes with its surroundings (Thomson and 
Newhall 1885). However, a dense drop of fluid with 
surface tension can avoid disruption and settle indef- 
initely at terminal velocity (as do raindrops, or water 
in oil). In this situation, the criterion for stability of a 
fluid droplet, which determines its terminal velocity and 
thence its size, is given by the Weber number We, where 
We = 2rp g v 2 /'-f, and 7 is the surface tension coefficient. 
Drops arc disrupted when We exceeds some critical value 
We* - 1 - 10 (Pruppacher and Klett 1997). 

Dense particle clumps in the nebula have no surface 
tension, but they do have self-gravity. By direct analogy 
with the droplet surface tension criterion, we define a 
gravitational Weber number Wee as the ratio of the ram 
pressure force per unit area to the self-gravitational force 
per unit area for a flattened disk, which is the initial 
stage of a strengthless spherical clump of diameter I and 
particle density p p = <&p g upon encountering a headwind 
w g (Thomson and Newhall 1885). For a disk of surface 
mass density a p — p p l 7 the gravitational force per unit 
area is Go 2 . Then 



We G 



C D p g w 2 C D w 2 



2Gc 



2G<$> 2 Pn l 2 



(2) 



where Cd is an effective drag coefficient for the clump, 
on the order of unity (see Appendix B). Wee can be 
written in other ways as well. The orbit frequency at a 
provides the useful relationship f2(a) = (GM a /a 3 ) 1/ ' 2 = 
(Gp*) 1 / 2 , or G = VL 2 /p\ where p* = M Q /a 3 ~ 3.8 x 
10" 8 g cm~ 3 at 2.5 AU 6 ; then 

C D p*w 2 g C D f3 2 p* a 2 



We G = 



2<$> 2 p g Vl 2 l 2 



Pg 



(3) 



6 Note p* is different from Safronov's (1991) Roche density pn 
3Mq /Atto? discussed earlier 



By analogy with the more familiar surface tension case, 
there will be some critical gravitational Weber number 
for stability Wee*, that is probably on the order of unity 
(Pruppacher and Klett 1997); its value must be con- 
strained by numerical experiments as described below. 
Then we require for stability against headwind disrup- 
tion that We G < We G or 

*>>^CWe M V>= {2a J£ /CD)1/2 . (4) 

These expressions determine the combination of clump 
diameter I and particle loading $ that will stabilize it 
against a headwind of magnitude fiv^- The headwind 
due to nebula radial pressure gradient will occur even if 
a clump is at the midplane and has zero settling veloc- 
ity; this gives the lower limit on ram pressure that the 
clump must be stable against. Neglecting vertical set- 
tling restricts potentially stable clumps to lying within 
about 0.01 gas scale heights of the midplane, where the 
headwind [3Vla is comparable to the vertical settling ve- 
locity. This restriction must be factored into statistical 
estimates of planetesimal production (see, e.g. Cuzzi et 
al 2007). 

For typical nebula parameters, we assume gas den- 
sity to be in the range p g = 10~ 10 g cm~ 3 (a nom- 
inal minimum mass value) to 10~ 9 g cm~ 3 (a more 
recent value supported by nebula evolution and chon- 
drulc formation; cf. Desch et al 2005), semimajor axis 
a = 3.8 x 10 13 cm (2.5 AU), pressure gradient/headwind 
parameter (3 <~ 10~ 3 (sec, eg., Nakagawaet al 1986, Cuzzi 
et al 1993), and solar mass M Q = 2 x 10 33 g. In this 
regime, $i > 1.5 - 5 x 10 6 (/3/10~ 3 ) km. Recall from 
section 2.3 and Appendix A that mass loading limits the 
maximum achievable mass loading to $ ~ 100. A clump 
satisfying the above constraint, with I = 1 — 5 x 10 4 km, 
would have the mass of a 10-100 km radius body of unit 
density - thus, these precursors can lead to quite size- 
able objects. This characteristic size range is intriguingly 
close to the roughly 50 km radius primordial building 
block size of Bottke et al (2005). The very existence of 
any preferred size for the primordial population, if true, 
is an intriguing result. 

3.4. Numerical model of clump evolution 

To obtain a sanity check on the concepts of section 
3.3, and to constrain the (unknown) value of We G for 
this problem, we developed a very simplified numerical 
model of a clump experiencing a steady nebula headwind 
from a more slowly orbiting, pressure-supported gas, us- 
ing Hill's approximation which transforms a cylindrical 
system into a cartesian system rotating at some mean 
rate fio - useful if the domain covers a sufficiently narrow 
radial and angular region that the curvature is negligible. 
The overall orbital motion of particles (the Keplerian ve- 
locity) is thus subtracted out, and the gas exhibits a 
small differential velocity because of its deviation from 
Keplerian. A particle clump under influence of the gas 
will slowly lag behind in the y (orbital) direction, and 
slowly drift to smaller x (radial) values, due to the gas 
headwind. The model is not intended to be a high-fidelity 
representation of a realistic nebula situation, but merely 
to check the basic Weber number model of section 3.3. 

The code evolves the velocity field of a two-phase sys- 
tem. The gas is perturbed by a clump of particles, which 
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are mutually attracted by gravitational forces and are 
themselves dragged by, and exchange momentum with, 
the gas. The Hill frame, with cartesian coordinates 
(x, y, z), represents a radially narrow region correspond- 
ing to some range of radius a, orbital angle 0, and ver- 
tical distance Z, respectively. Periodic boundary condi- 
tions are used in all three dimensions. The instantaneous 
Navier-Stokes (Eulerian) equations describing the con- 
servation of mass and momentum for the incompressible 
gas are expressed in the Hill frame as, 



V-U = 0, (5) 



and 



^ + ( u.v)u= VP+Vpg ' o6 + ,vnj 

dt p g 

-2n x U + i drag + i sp (6) 

where U is fluid velocity, p g is gas mass density, v is 
gas viscosity, P is local pressure (which fluctuates along 
with gas velocity and vorticity variations), \7P 9lob = 
(2p g WgQo,0,0) represents the (constant) global nebula 
radial pressure gradient, and fi = (0,0, £7q) is the Ke- 
pler frequency vector. This setup causes the gas to 
flow towards the clump in the y direction at a uniform 
speed w g . We are currently neglecting the Keplerian 
radial shear term 2qQ 2 xe x , where f2(a) = flo(a/a )^ q 
and Oo = f^(ao), in the gas and particle equations 
to get a better comparison with our analytical stabil- 
ity model. The shear rate associated with this term 
(dw / dx)xep = whereas a typical (fluctuating) 

local shear or strain due only to turbulent motions is 
(dw/dx)turb ~ > O . Because the turbulent shear 
is so widely varying, and likely to be larger than the 
Keplerian shear, in this study we neglect both of these 
complications, though they are both suitable avenues for 
future work. 

In our code, particles are followed in the Lagrangian 
sense. The term f dra 9 represents the force per unit mass 
imparted to the gas by the particles. It has the general 
form 



^drag 



Pp 
Pgte 



(V-U) 



(7) 



where V is the mean weighted particle velocity at a grid 
point, t s is the particle gas drag stopping time, and p p 
is the particle mass density. To obtain V, a weighted 
sum is carried out over particles in the eight cell volumes 
adjacent to each fluid grid point. Weighting functions 
are used which vary inversely with the distance of the 
particle from the grid point (Squires 1990, his section 
5.1). Because of the periodic boundary conditions, the 
wake of a clump can impinge artificially on the clump 
from the upwind direction. To avoid this, we include the 
term f sp in equation (6) as a "sponge force" per unit mass 
to restore the gas velocities downstream of the clump at 
the outflow boundary plane to their initial values U ; this 
ensures constant inflow values at the upstream boundary 
plane. 

The sponge force has the form 



f sp = — (U-U ), 



(8) 



where r sp is some time scale. The function 1 / T sp is mod- 
eled as a sigmoid in the y direction very near the outflow 
boundary. The parameters of the sponge were deter- 
mined through test runs with the goals of minimizing 
the sponge's spatial extent and maximizing its effective- 
ness without introducing numerical instabilities. Once 
determined, the same values were used for all produc- 
tion runs. 

The (Lagrangian) equation describing the motion of 
particle i subject to the forces of gas drag and mutual 
gravity is (again neglecting the Keplerian shear term) 

dVi 



dt 

pdrag 



2fi x V, 



pdrag 



ff 



(9) 



The term ^ 9 describing the drag force per unit mass 
on particle i by the gas takes the form 



pdrag 



i(U(X t (t))-V,(t)), 



(10) 



where U(Xj(£)) is the gas velocity interpolated to the 
particle's position Xj(f) at time t. 

Finally, the mutual gravity force per unit mass on each 
particle is 



cgrav 



N p 



X«(t) 



(|(Xy(i)|+e)' 



(11) 



where G is the gravitational constant, and N p is the num- 
ber of particles. Xj 3 is the unit vector associated with 
the distance Xy (t) = Xj (t) — Xj (t) between particles j 
and i, and e is a constant that softens the force at small 
separations to prevent numerical singularities, chosen to 
be comparable to a grid cell in extent. 

Eqs. (6) and (9) are solved using psuedo-spectral meth- 
ods commonly used to solve Navier-Stokes equations for 
a turbulent fluid (Canuto et al 1987). Periodic boundary 
conditions are assumed and the number of nodes used in 
each direction is generally 128, 256, 128 with a spacing 
of 27r/64. A Fast Fourier Transform (FFT) algorithm is 
used to evaluate the dynamical variables U at the com- 
putational nodes. A second-order Runge-Kutta scheme 
is used to time-advance the gas and particle velocities. 
A third-order Taylor series interpolation scheme is used 
to determine gas velocities at the particle positions from 
values at the eight nearest neighbor nodes. The mutual 
gravity calculation is done in a brute force fashion with 
a N 2 algorithm to evaluate the force contributions from 
all particles. The code is written in Fortran 77 and is 
parallelized using OpenMP directives. The FFT calcula- 
tions are done simultaneously on the planes perpendicu- 
lar to the y and z axes, and all loops involving particle 
indices are parallelized. Typically, 3000 "supcrparticles" 
are used, but some runs used more than 20000. Each "su- 
perparticle" represents the dynamical effects, and typical 
response, of a large number of actual chondrulcs. 

Initially the particles are arranged in a uniform spher- 
ical clump. The initial velocities for the gas and the par- 
ticles inside and outside the clump are determined from 
the expressions in Nakagawa et al. (1986): 

p p / 2L>0 D 2 p g + p p 



U = Wg 



V = -Wg 



p p \D 2 + n 2 ' D 2 + n 2 

2DQ D 2 



Pp 



Pg 



Pg + p p \D 2 + n 2 ' D 2 + n 2 



,o 



Am 



(13) 
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Fig. 3. — Snapshots from evolutionary models of dense clumps 
experiencing a nebula headwind without (left column) and with 
(right column) self-gravity. Top row: projected into a vertical plane 
(z—8); bottom row: projected onto the a—z plane (a— 8 projections 
are similar), with the gas inflow into the page. In the top left 
panel the open circle shows the initial size and position of the 
clump. After less than an orbit, ram pressure disrupts the gravity- 
free clump (left) but the gravitating clump is stable (right) and 
continuing to shrink inexorably, even while being eroded. Movies 
showing this evolution arc available online (see Appendix C). 

where D = (p g + p p )/ p g t s . The initial velocity of the gas 
outside the clump is found by setting p p — > in equation 
(12): Uo = (0, —w g , 0); inside the clump the specific case 
value of $ is used in equations (12) and (13) to initialize 
the gas and particle velocities. 

Runs with single particles of different t s were made 
to verify that the Nakagawa et al. (1986) initial condi- 
tions are steady state solutions for the particle equations. 
The wallclock time to evolve one integration step is 0.7 
sec using 256 Intel Itanium 2 processors running at 1.5 
gighertz. The wallclock time to evolve 1 orbital period 
is 121 hours. 

3.5. Results of the numerical model 

Figure 3 shows how, as predicted by the simple theory 
of section 3.3, certain combinations of clump mass load- 
ing <I> and dimension I remain stable against ram pressure 
disruption for nebula headwinds produced by a pressure 
gradient characterized by (3. We is dimensionless, and 
the second expression in eqn. (2) can be used to obtain 
the critical value of Wee which separates stable from un- 
stable configurations of the clump, Wee = We^, directly 
from measured values in code units (see Appendix B). 
We expect that We£j will be in the range 1-10. The com- 
bination of parameters in the stable case shown (right 
panels) gives We^j ~ 1; if self gravity is turned off, the 
clump is disrupted in tdis ~ 1/27T orbits, whereas it will 
sediment into itself on a timescale t se d ~ 8 times as long 
(see Appendix B). In the stable cases, a dense core is seen 
to continually shrink and become denser throughout the 
run, even as material is shed from the periphery of the 
clump, such that the value of <&Z for the core continues 
to increase (see Appendix C, figure 7). 

The numerical clump models exhibit considerable ero- 
sion over the duration of the runs from a viscously stirred 



surface layer, and this gradual but inexorable mass loss 
might lead to their disruption if arbitrarily long runs were 
practical at this time. However, as discussed in Appen- 
dices B and C, this large amount of erosion is a gross over- 
estimate of what would happen in the nebula, because 
numerical viscosity (and consequently shear around the 
periphery of the clump) plays a far larger role in the cur- 
rent numerical models that it ever would in the actual 
nebula. That is, the viscously stirred and eroded surface 
layer in the numerical clumps contains orders of mag- 
nitude more mass than would be the case for a nebula 
clump of the size and density required to survive ram 
pressure disruption (see Appendix B). In Appendix C, 
we describe the movies from which figure 3 is taken and 
also show how the protected inner regions of the clumps 
are behaving exactly as predicted, both in stable and un- 
stable regimes. For more realistic, larger clump Re c in 
the nebula, a far larger fraction of the clump mass will 
show this behavior rather than being artificially eroded. 

The primary purpose of the numerical models is merely 
to provide an independent sanity check on the physics 
of our Weber number model and obtain some idea of 
We^, which gives us an order-of-magniturc estimate for 
the product <&l which nebula clumps require to survive 
headwinds of a particular f3p g . Future coding improve- 
ments are needed to allow larger (higher Re c ) clumps 
which would more closely approach nebula conditions in 
terms of their balance between pressure forces and vis- 
cous forces (see Appendices B and C); these improve- 
ments will probably include an implicit time advance for 
the drag terms and perhaps also a tree or particle-in- 
mesh code for the particles. 

4. DISCUSSION AND SUMMARY: 

We have shown how self-gravity can stabilize dense 
clumps of mm-sized particles, which form naturally in 
3D turbulence, against disruptive gas ram pressure on 
timcscales which are sufficiently long for their constituent 
mineral grains to sediment towards their mutual cen- 
ters and form physically cohesive "sandpilcs" of order 
10-100km in size. The essence of the result is a criti- 
cal "gravitational Weber number" on the order of unity, 
in which self-gravity plays the role of surface tension in 
more familiar situations such as raindrops. Characteris- 
tic mass densities and lengthscales are determined which 
meet this requirement. We show numerical results which 
are in general agreement with the predictions of the sim- 
ple theory, for isolated, spherical initial clumps. 

The scenario we have sketched out leads from aerody- 
namically size-sorted nebula particles, having the prop- 
erties of chondrulcs, to sizeable planetesimals formed en- 
tirely from such particles which contain a snapshot or 
grab-sample of the local particle mixture. Turbulent con- 
centration first produces the dense zones of size-sorted 
particles. The more common, less dense of these re- 
gions may provide typical chondrule melting environ- 
ments (Cuzzi and Alexander 2006). Some of the less 
common, very dense zones, which we have shown else- 
where can achieve mass densities 100 times larger than 
that of the gas, have the potential to become planetes- 
imals depending on their lengthscales, nebula location, 
and local vorticity. It is intriguing that our character- 
istic stabilized clump masses are not far from the mass 
inferred by Bottke et al (2005) to represent a typical pri- 
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mordial object in the pre-dispersal, prc-erosional asteroid 
belt. 

Ultimately, the sandpiles resulting from completed sed- 
imentation will become compacted further by inevitable 
collisions with other sandpiles, leading to today's fairly 
dense asteroids, while retaining a physical and chemical 
memory of their parent particle clumps. The mechanism 
is easily extended to the unmelted aggregate particles of 
the outer solar system, where, however, if chondrules are 
absent, the size-sorting fingerprints of the process (Cuzzi 
et al 2001) might be less evident. 

The conclusions of this paper differ from the sugges- 
tions in Cuzzi et al (2001), who focussed on the possible 
role of ultra-dense clumps with size comparable to a Kol- 
mogorov scale (0.1 -1 km). Since that time, Hogan and 
Cuzzi (2007) found that particle mass loading saturates 
the value of $ = p v j ' p g at a value of about 100, redirect- 
ing our attention to larger clumps in the 10 3 — 10 4 km size 
range. Clumps and fluid structures of these large sizes 
are much more accessible with numerical codes, both of 
the standard direct simulation type and the cascade type 
(Appendix A). 

In this scenario, primitive bodies may not primarily 
represent spatial, but perhaps temporal, samples of the 
particulate contents of the nebula as its chemical, physi- 
cal, and isotopic properties evolve over several Myr. The 
temporal, as well as the spatial, variation in the physi- 
cal, chemical and isotopic properties of the concentrated 
particles, which are being continually altered by thermal 
events and mineralogical alteration in the nebula gas, can 
then help account for class-to-class variations between 
the chondrite groups. An implication of this drawn-out, 
inefficient process is that younger chondrite types should 
contain evidence of "leftovers" or "refugees" from earlier 
times, which escaped primary accretion. Several aspects 
of chondrite makeup are compatible with this implica- 
tion. Of course, it is well known that ancient CAI min- 
erals and less-ancient, less refractory Amoeboid Olivine 
Aggregates are found alongside much younger objects in 
the same chondrites (eg. Scott and Krot 2005). Also, the 
oldest chondrites (CVs) contain primarily type I chon- 
drules (which have an age observationally indistinguish- 
able to that of CAIs) and no (more oxidized, heavier 
O-isotope) type II chondrules, while the younger ordi- 
nary chondrites contain a mixture of type I and type II 
chondrules, some with hints of age variation even within 
a given chondrite (Kita et al 2005). The oldest and 
first-formed objects are all likely to have melted from 
their abundant 26 Al (Kleine et al 2005, 2006; Hevey and 
Sanders 2006), producing differentiated achondrites and 
metallic objects, so it is no longer possible to dissect their 
primordial components. 

Of course, the real world is more complicated. Some 
complications are meteoritic: as only one example, CO 
chondrites and ordinary chondrites appear to be about 
the same age (younger than CV chondrites), but have 
very different chemical and isotopic properties (Kita et 
al 2005), which would, in the context of this scenario, 
argue for some degree of spatial (radial?) variation in 
the makeup of nebula particulates at that time anyway. 

On the theoretical side, of course, self-consistent nu- 
merical models which follow clump evolution in realis- 



tic turbulence with headwinds and vertical settling need 
to be pursued, to check these preliminary assessments 
which make simplistic assumptions about the local head- 
wind and assume isolated clumps. Several recent stud- 
ies have found that transient "pressure ridges" arise in 
the gas in association with spiral density waves or vor- 
tical structures (Haghighipour and Boss 2003, Rice et 
al 2004, Johansen et al 2007); in such regions, the local 
radial pressure gradient and headwind diminish (other 
regions will have unusually large headwinds). Moreover, 
clumps do not exist in isolation, but more realistically as 
a dense core within a larger, less dense envelope which 
weakens the headwind felt by the dense clump core rel- 
ative to the nebula average we characterize by (3. Both 
of these conditions would allow some clumps to survive 
with smaller $Z. Other complications include interac- 
tions between different strengthless objects, leading to 
mergers or disruptions. 

Our stability criteria are most likely to be satisfied in 
a region fairly close to the nebula midplane (~ O.OliJ, 
see section 3.3), because at higher altitudes the verti- 
cal component of solar gravity leads to vertical settling 
of dense clumps at higher terminal velocities than are 
easily stabilized against. However, this vertical settling, 
even if it leads quickly to disruption of individual clumps, 
enhances the downward transport rate of small parti- 
cles and increases the solid/gas ratio closer to the mid- 
plane from that generally predicted by simple ID diffu- 
sion models (eg., Bosse et al 2006). 

All these effects have implications for the occurrence 
frequency of clumps of the requisite <frl for stability. 
Quantitative statistical determinations of the volume 
fraction of stable clumps which can become planctcsi- 
mals at any given time will use, for instance, approaches 
such as the cascade models described in Appendix A. 
Simplified, preliminary work exploring these factors in- 
dicates that, while the volume fraction of stable clumps 
is low at any given time (so the process is clearly not 
an efficient one), accretion rates are roughly an Earth 
mass per Myr in the asteroid belt region (Cuzzi et al 
2007). However, considerable work remains to establish 
the statistical formation rate of appropriate clumps ca- 
pable of following this evolutionary path and to study the 
evolution of clumps in realistic turbulence, with variable 
headwind velocity and simultaneous vertical settling. 
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Fig. 4. — Conditioning curves for multipliers p(m); the parameter b depends on the local mass loading larger b indicates a narrower 
PDF p(m). Open symbols: for particle concentration C; filled symbols: for cnstrophy (S = ui 2 ). The small insets show p(m) at several 
values of <t>. From Hogan and Cuzzi (2007). 




Fig. 5. — Joint (binned) probability <&u) 2 P(<fr,uj 2 ) from the conditioned cascade model, for a 24-level case with initial Pp/Ps = 10~ 2 . 
The curves of figure 1 are the S— averaged equivalent of this plot (and show C instead of <3?), but do not include the effects of mass loading. 
Note the flattening or saturation at <E> ~ 100, reflecting the near-vertical asymptote of the conditioned-fe curve of figure 4 at <I> ~ 100. A 
24-level cascade is in the plausible range for clumps of size ~ 10 4 km in a nebula with a ~ 10 -3 — 10~ 4 . In this figure, the enstrophy S is 
normalized to its mean value for the binning lengthscale (from Hogan and Cuzzi 2007). 

APPENDIX 

APPENDIX A: MASS LOADING AND THE CASCADE MODEL 

This aspect of our work incorporates two related and important concepts: that of intermittency, and that of a 
statistical cascade process. For instance, it is widely known that dissipation of turbulent kinetic energy occurs at 
the Kolmogorov scale; it is less widely known that the spatial distribution of dissipation, like that of the particle 
concentration factor C, is highly intermittent. Intermittent quantities are spatially and temporally unpredictable, 
and fluctuate increasingly with increasing Reynolds number Re. However, the statistical properties of intermittent 
quantities like C (their probability distribution functions or PDFs) are well determined on any lengthscale. In what 
follows, it will be necessary to distinguish between two closely related quantities: the concentration C and the mass 
loading $ = p p /p g = C~p^/p g . This is because the emergence of particle feedback on the gas due to mass loading ($), 
depends both on C and on the initial value of p^. It is ultimately C that is evolved in the cascade, but how it evolves 
depends on $. So both quantities will be alluded to in parallel. 

In its inertial range, which is extensive at high Re, turbulence is a scale-free process that is often referred to as a 
cascade and inspired the approach of cascade modeling (Meneveau and Sreenivasan 1991). The physics of transport 
of kinetic energy, vorticity, and dissipation from their sources at large scales is independent of scale until viscous 
processes enter at the smallest (Kolmogorov) scale. Scale-independence and Re— independence are connected because 
Re determines the depth of the inertial range: Re 3 ' 4 = L/r) (Tennekes and Lumley 1972). Cascade models simplify 
this complex, nonlinear, 3D scale-free process into a set of partition rules which are independent of scale, or of level 
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in the cascade. Larger Re means a larger number of eddy spatial bifurcations (or levels) between L and r), and 
stronger fluctuations in intermittent properties (Meneveau and Sreenivasan 1991). Cascade models don't preserve 
all the information of full 3D models (such as the tubelike spatial structures characterizing vorticity) but they do 
reproduce the statistical properties (the PDFs), which are of primary importance for our problem. The successful use 
of such models to describe the intermittency of dissipation, and other quantities, in turbulence (Juneja et al 1994, 
Sreenivasan and Stolovitsky 1995) led us to pursue a cascade model for turbulent concentration of particles. 

A cascade model consists of a set of "multipliers" - fractions which define how a quantity of interest is unequally 
partitioned from a fluid parcel at some level, into its equal-sized subdivisions at the next level. The PDF of particle 
concentration factor C is determined over all the numerous end-level outcomes of applying a multiplier m (chosen 
from its own PDF p(m)) to each volume clement as it bifurcates at each level of the cascade. It has been observed 
that p(m) is independent of level in the turbulent inertial range (Sreenivasan and Stolovitsky 1995). Hence, multiplier 
PDFs determined from a direct numerical simulation (DNS) at low Re (with a relatively small number of levels), 
merely applied repeatedly over additional levels, can predict the properties at higher Re (if nothing changes in the 
physics). Our cascade model (Hogan and Cuzzi 2007) uses a common form for p(m): p(m) oc m b_1 (l — to) 6 ^ 1 , where 
the parameter b determines the width of p(m) (Sreenivasan and Stolovitsky 1995). Both the concentration C and the 
vorticity w are described by their own ^-distributions. Multiplier PDFs with small b(< 1) are broad, and in them 
multipliers much different from m = 0.5 occur with higher probability, producing a grainier, more intermittent spatial 
distribution. PDFs with large b(^> 1) are narrow and centered on m = 0.5; because each "eddy bifurcation" then 
involves a nearly 50-50 partitioning, the resulting spatial distribution is nearly uniform and subsequent growth of $ 
becomes more difficult. 

The cascade model can then be used to address the issue of how particle mass loading can affect the cascade - that 
is, change the physics as the cascade progresses. This dependence is known as conditioning of a cascade. We found 
that the process can be represented quite well as two separate one-phase conditioned cascades for C and u> in which 
multiplier distributions p(m) for both C and lu depend only on the local particle mass loading $. We used our full 3D 
numerical simulations to establish how b for both these properties depends on local particle or fluid properties. 

Figure 4 shows conditioning curves for both C and u> (the latter expressed as enstrophy S = uo 2 ), as extracted from 
multiple 3D DNS simulations at various Re as large as 2000 (Hogan and Cuzzi 2007). Note how the multiplier PDFs 
extracted from the DNS results (shown using insets) get narrower (have larger b) as mass loading increases, choking 
off intermittency. In these runs, ~p^/p g = 1, so C = $. These results establish the upper limit <f> ~ 100 which can 
be obtained by turbulent concentration, regardless of the initial value of ~fo/p g - 111 the nebula, with a smaller initial 
J^/pg, higher Re or deeper cascades (and larger ensuing C values) are needed to reach this saturation point. 

To determine the joint PDF of concentration C and vorticity u>, we require not only the conditioned multipliers for 
u 2 and C (Figure 4), but also their spatial correlation. Hogan and Cuzzi (2007) found, on average, a 70-30 preference 
for anticorrelation at each partitioning, consistent with previous observations that particle concentration zones avoid 
zones of high fluid vorticity (Squires and Eaton 1990, 1991; Eaton and Fessler 1994; Ahmed and Elghobashi 2000, 
2001). Including this partitioning asymmetry factor as a weighted coin gives us a two-phase (particle-gas), conditioned 
cascade model which shows very good agreement with the joint PDF P{<&,uj 2 ) directly determined from our full 3D 
mass loaded simulations (see Hogan and Cuzzi 2007). However, to match the full 3D DNS simulations at Re = 2000, 
the cascade model needs only about 15 levels, taking about 10 cpu-hours (for 1024 realizations) compared to over 
90000 cpu hours to converge our full 3D simulation. 

Using these conditioned multipliers and asymmetry factor, entire PDFs of particle concentration (mass loading) and 
vorticity can be generated for arbitrary numbers of levels (arbitrary Re). For example, the 24 level model of Figure 5 
(from Hogan and Cuzzi 2007) is of direct relevance for 10 3 — 10 km scale dense zones of interest for the nebula. Figure 
5 illustrates the saturation of the particle concentration at <I> = 100 and the fact that high <I> occurs preferentially at 
low enstrophy or vorticity. The upper left quadrant is of most interest for determining the numbers of clumps which 
can survive to become planetesimals (see Cuzzi et al 2007 for more details). 

APPENDIX B: SCALING BETWEEN NUMERICAL CODE AND NEBULA 

While the primary result of this paper - the existence of some stability regime determined by a critical gravitational 
Weber number We^ - can be expressed in a nondimensional way, some aspects of the results (specifically the numerical 
results and the inferred value of $1 for the nebula) require us to pursue the relationship between nebula parameters 
and code units more deeply. Values in code units are denoted by primes below. The fundamental quantities in the 
problem are (1) the gas and particle densities p g and p p = Qp g ; (2) the clump diameter I; (3) the local orbit frequency 
O(a) where a is the distance from the sun; (4) the local radial pressure gradient dP/da = dP/dx, where x is in the Hill 
frame; (5) the particle stopping time t s ; (6) the velocity of material in the Hill frame W, that is, relative to Keplerian, 
where W represents either the particle or gas velocity; and (7) the gravitational constant G. The principal forces are 
the coriolis force due to the rotating frame — 2f2o x W, the pressure force (— 1/ p g )dP/dx, and the gas drag forces that 
couple the gas and particles (section 3.3). The code unit of length is radians (the domain is 2ir radians on the short 
dimension (x or z) and the grid cell size is therefore A = 2tt/N where N is the number of grid points along that axis). 
The code time unit (ctu) is where O' = 1 is the code value of rotation frequency; thus the orbit period is 2n in 
code time units or ctu. We treat as equivalent a = R = x. 

As noted in section 3.3, 



n(a) = (GM /a 3 ) 1/2 = (Gp*) 1/2 , or G = Q 2 /p* 
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where p* = A'/ Q /a 3 ~ 3.8 x 10~ 8 g cm~ 3 at 2.5 AU (note this is different from Safronov's (1991) Roche density 
Pr = SMq/Atto 3 ). To normalize mass and mass density (and ultimately establish the code value of the gravitational 
constant G') we assume a nebula gas mass density p ~ 10~ 10 g cm~ 3 (minimum mass nebula; Cuzzi et al 1993) at 
2.5 AU; then code quantities p' g = p g /p = 1 and p' p = p p /p — ^p' g — < f ) , and p'* = p*/p - Also, in code units 
O' = £!/Oo = 1 where ft = 0(a o ) is the actual orbit frequency at the region of interest. Then the gravitational 
constant in code units, G", can simply be written as G = Q' 2 / p'* = 1/ p'* ■ Finally, the pressure force is written as 

1 dP 

— = 2/M2 2 = 2/3aGp* = 2w q fl 1 

p g ax 

where j3 « 10~ 3 is typical (Nakagawa et al 1986; Cuzzi et al 1993), G' — p /p* — 3 x 10~ 3 . Then the dimensionless 
parameters $ and /' are set to satisfy various constraints of the code (resolvable clump, adequate particle statistics, 
reasonable timescale, etc) and (3' is allowed to vary, to determine empirically the critical value of We^. as the largest 
value of Wee that remains stable. That is, we want /' to be large enough such that the clump is well resolved, but not 
so large that it fills the entire cross section of the computational box. Typically /' is only 10-20 grid cells so far, and 
as noted below, this exaggerates the viscous stresses and surficial erosion of material. We want $ to be large enough 
that the clump is much denser than the gas, and we want a large enough number of particles to act like a continuum. 
The value of t s should ensure that the clump sedimentation time t se d is larger than its ram pressure disruption time 
tdis (see below) , to provide a real test that self-gravity is preserving the clump rather than simple collapse. 

Important timescales in the numerical model 

The dynamical collapse time of the clump under its own self-gravity, and in the absence of gas pressure, is to = 
Tr(Gp p )~ 1 / 2 , or in code units t' G — 7r(G"$) -1 / 2 . The mass loading which produces a clump dynamical time comparable 
to the orbit time is then obtained using 7r(G' c E > )~ 1 / 2 = 27r/0' = 2n or $ = 1/4G" ~ 100, again in fair agreement 
with prior expectations. However, gas drag prevents clump collapse on this timescale, as described in section 3.1, and 
actual shrinkage takes a time t se d (eqn. 1). The number of orbits over which we need to follow the clump to ensure it 
really survives is roughly 

t sed n cVt Hfl 2 fl 2 H 1 _ p* H 1 

27r 8TrG§rp s 8TrG§rp s Gp s r 87r$ p s r 87r$' 

where we have used c = HVl. This expression can be assessed using real quantities, and is 30-300 orbits for particle 
radius 300/x and mass loading $ = 1000 — 100. This behavior (t se d S> 2ir/n') can be controlled in the code, once 
the other parameters are established as above, by adjusting the particle stopping time t s . Because of long code run 
times, at present we are limited to stipulating only that the sedimentation time significantly exceeds the nominal 
(self-gravity- free) disruption time; that is, t se d ^> tdi S , where 

t dis « (l/w)y/2*/C D . 

This expression for the disruption time tdis is derived from the ram pressure force (sect 2.3) by setting the distance 
traveled in time t by the windward half of the clump, under acceleration by the ram pressure force, equal to the clump 
diameter (assuming the leeward half of the clump doesn't get accelerated), and solving for t. The value of Cd for our 
clumps is on the order of unity, as we verify below. 

Numerical viscosity and the clump Reynolds number: 

The code calls for some input viscosity u'; with it, the Reynolds number of the clump having diameter I' radians, 
in a headwind of speed w' g (in code units of radians per ctu), is Rc' c = l'w'/2v' where v' = 0.1 is the defined code 
viscosity (radians 2 /ctu). We assume an inertial range expression within the wake of the clump to determine the wake's 
Kolmogorov scale: rj'/l' = Rc' c ~ 3 / 4 , and thus rf = [2v'l ll l' i jw'gf'^ radians. Our highest resolution runs to date had 
a gas relative velocity w' g = 38 radians/ctu and used a clump size V = 2 radians. The nominal Kolmogorov scale 

associated with this flow is rf = (2(0.1)(2) 1 / 3 /38) 3 / 4 ~ 0.02 radians (compared to the grid cell size of 0.05-0.1 radians); 
the wake turbulence is thus under-resolved and the code's true viscosity is numerical: v' n ~ w' g A (at the boundary of 
the clump). 

While we are not concerned with the fine-scale details of the wake, this marginal resolution introduces some caveats. 
The ratio of ram pressure force to viscous force (both per unit area) is, in general, (p g w 2 /2) / (p g vdw g / 'dr) ~ w g ljAv = 
Re c /2 where Re c is the Reynolds number of a clump, and we have approximated dw g /dr ~ w g /(l/2). Because we are 
dominated by numerical viscosity, v' n ~ w' g A, thus in code units Re^. ~ w g l'/2w' g A = l'/2A « 5 in many cases so far. 
For comparison, the Reynolds number for nebula clumps of interesting sizes is > 10 6 . This means that viscous stresses 
around the periphery of the clump, and the fractional depth of the viscous boundary layer, are grossly exaggerated 
in the numerical model relative to the actual nebula regime of interest. The anomalously large role of (numerical) 
viscosity leads to anomalously large "erosion" from the surface of our numerical clumps, and degrades our results in 
the sense that potentially stable clumps might appear less stable, due to the erosive mass loss experienced over their 
sedimentation time. Going to higher resolution cases (Re^, ~ 10) significantly ameliorated this effect but the problem 
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Fig. 6. — Determination of the effective drag coefficient for our clumps, from plotting estimated disruption times t^ is vs. the quantity 

(//w; g )(2<l>) 1 / 2 ; the slope is then C^, 1 ^ 2 . This determination is slightly qualitative as to when we judge the clump to be disrupted; we 
decided to use the time at which vortex-driven voids first appear in a a — z planar cross-section through the nominal center of the clump. 

has certainly not been entirely removed. Another way to characterize the degree of artificial erosion is to estimate the 
physical Kolmogorov scale for the actual nebula/clump situation: for v = 10 cm sec" - 2 , I = 10 4 km, and w g =3800 
cm/s, we find rj ~ 100m, or 10~ 5 of the size of a 10 km clump, compared to a fractional size of perhaps 10 _1 or so 
in the code, as determined by numerical viscosity. Clearly, surface erosion is a much smaller effect in the nebula case 
than in our crude models. While even at this resolution apparently stable configurations can be found, it would be 
desirable for future studies to improve on this situation. 

Clump drag coefficient Cd ■' The drag coefficient of the clump is Reynolds- number-dependent and must also be scaled 
to nebula values to obtain an estimate of $Z. Weidenschilling (1977) gave expressions for the drag force per unit area 
(giving the pressure force per unit area), experienced by solid spheres of different Re c in the Stokes drag regime of 
interest here, as CoPgWg/2. For 1 < Re c < 800, Cd — 24/Re°' 6 , or Cd ~ 8 for Re c ~ 5, as in our code. The large 
drag coefficient and ram pressure is partly due to low-pressure zones set up immediately behind particles in this range 
of Re c . By comparison, for a solid sphere with Re c > 800 as under nebula conditions, Cd = 0.44. However, neither of 
these values might be entirely applicable to our clumps, which are not rigid spheres. In practice, we have determined 
Cd empirically from observed values of tdis with self gravity turned off. Recalling that tdis = (I /w g )(2$ / 'Cd) 1 ' 2 ', we 
simply plot our estimate of tdis, the point at which sizeable vortex-driven internal voids appear in the clump, against 

the combined parameter (Z/u> g )(2$) 1 / 2 , giving C^, 1 ^ 2 as the slope of the best fit straight line (figure 6). It seems that 
Cd ~ 1 is a good assumption, but this calculation should be redone with higher resolution (lower numerical viscosity) 
codes. At nebula Re c , Cd could reach its theoretical high-i?e value of 0.4 and this is what we assume in deriving <I>L 
Finally then, the several stable cases we have run (eg. figure 3) can be characterized by w' g = 38 rad/ctu, p'* — 380 

(which corresponds to a minimum mass nebula local density p g = 10~ 10 at distance 2.5 AU), <E> = 10 3 , and I' « 1 
radian, and while we have not as yet accurately determined the transition value of Wee, a value Wee ~ 1 is clearly 
stable. Solving equation (4) of section 3.3 for $1 under nebula conditions, assuming We^j = 1, taking Cd = 0.44 for a 
high-Re c situation and assuming a nominal j3 — 10 -3 , we find $1 > 1.5-5 x 10 6 , independent of p g and neglecting a 
questionable coefficient of order unity (because of the contribution of viscous stresses in our low-Re c code). Foreseeable 
refinements (specifically, lower numerical viscosity) may increase We^ ( the surface tension analogue gives We^ ~ 10) 
and thus decrease the value of $L 

APPENDIX C: TIME EVOLUTIONS AND MOVIES 

Quicktime movies showing our time evolutions are available online at 
http: / /spacescien ce.arc.nasa.gov/users/cuzzi/| 

with filenames as given below. After publication, the movie files will also be available on the ApJ web site as 
mpgs. Each shows a dense clump in the Hill frame, rotating at local Keplerian velocity, with more slowly-orbiting gas 
impinging on the clump from its leading side. The two left-hand panes show the view from along the radial axis (in 
the Z — 8 plane; top) and along the vertical axis (in the a — 8 plane; bottom). The large right-hand pane shows the 
view from the orbitally leading direction (in the a — Z plane). The wake formed by the clump is easily visible in the 
two left-hand panes. Cases hill26.gif and hill29a.gif are sampled at a similar time near their end-point in figure 3. File 
hill26.gif shows how a clump with no self-gravity is disrupted in about 1.5 code time units for the initial conditions 
chosen (the orbit period is 2tt code time units), in agreement with our simple estimates of tdis- Files hill29a.gif and 
hill31a.gif are intended to illustrate the same (stable) Weber number case (Wee ~ 1) but we vary two constituent 
parameters of Wee (we decrease both the gas velocity and mass loading by a factor of three) to illustrate similarity; 
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Fig. 7. — Properties of the dense central regions in three simulations, as a function of time. Cases are shown in each panel as: solid 
line (stable, self-gravitating case 31a); dashed line (stable, self-gravitating case 29a); and dotted line (unstable, non-gravitating case 26). 
Not all cases were run to the same stopping point due to resource constraints. Top: cumulative total mass retained by the clump; middle 
panel: average mass density of the central mass fraction /; bottom panel: effective radius containing mass fraction /. These results 
show that, away from the artificially perturbed perimeter regions, W eQ-stablc clumps are behaving as expected by the simple model, and 
Wee-unstable clumps do not show this behavior. 

because of the specific parameters chosen, hill31a.gif runs for a longer time, but is suffering considerable erosion 
towards the end of the run. Nevertheless, as discussed below, it retains a dense core that continues to contract. As 
discussed in Appendix B, more realistic cases, with higher numerical resolution, would have smaller numerical viscosity 
and incur less erosion. Case hill29a.gif, while not run as long as hill31a.gif, appears solidly stable with a dense core 
that is continuing to shrink at the end of the run. Note some interesting oscillatory behavior, such as hinted at in our 
ID compressible runs (figure 2). 

Behavior of central regions of dense clumps 

The simulations which created the movies referred to above can be used to assess the behavior of the central, dense 
regions of each clump in a quantitative manner. Figure 7 has three panels, showing the time variation of some 
properties of the particles which lie within the central regions of the clump (its densest fraction, as sampled out to 
some cumulative mass threshold). In the top panel, we show how much mass is being eroded from the clump overall 
(particles moving faster than some comoving velocity threshold are assumed to have left the clump). All three clumps 
are losing mass, naturally, but the non-gravitating clump 26 is losing it much more rapidly than the gravitating clumps 
29a and 31a. In the central panel we show the mean mass density as measured over some fraction / of the clump 
mass lying at the highest densities; in clumps 26 and 29a, we chose / = 0.5, and in clump 31a, we chose / = 0.25. In 
the nongravitating clump, the core density never increases and it quickly blows up, while the central regions of the 
gravitating clumps (removed from the artificial mass erosion at their peripheries) inexorably get denser at a steady 
and perhaps even slightly increasing rate (an increasing rate is predicted by the simple model). The lower panel shows 
the effective average radius of the region being sampled. The non-gravitating clump never shrinks at all, but both 
gravitating clumps do. Moreover, as the gravitating clumps shrink, their density increases faster than their linear 
dimension decreases, so the product $>l increases - causing them to become more stable as time goes on. This is an 
argument for stability even in the presence of erosion. As expected, the density growth timescale (the sedimentation 
time t se d of equation 1) is three times smaller for the denser clump of case 29a (middle panel). Meanwhile, the apparent 
agreement of the variation of the fractional radius containing half of the mass seen between cases 29a and 31a in the 
lower panel is fortuitous, because the total masses of the clumps are changing (and eroding) at different rates and from 
different depths in the two cases. Nevertheless, the cores are both shrinking and becoming denser in the two stable 
cases. The amount of erosion that is incurred from the margins of each clump is related to (perhaps proportional to) 
the thickness of the viscous boundary layer relative to the radius of the clump, which can be expressed in terms of the 
Reynolds number. As we have argued above in this Appendix, actual nebula clumps capable of remaining stable would 
have Reynolds numbers five orders of magnitude larger - impossible to study with numerical models. The relevance of 
all viscous boundary layer effects, including the erosion that inflicts our current numerical runs, will decrease by that 
factor in the real situation. 
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